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Preface 

In  late  1976,  a  study  to  produce  a  wave  climate  for  U.  S.  coastal 
waters  was  initiated  at  the  U.  S.  Army  Engineer  Waterways  Experiment 
Station  (WES).  This  Wave  Information  Study  (WIS)  was  authorized  by  the 
Office,  Chief  of  Engineers,  U.  S.  Army,  as  a  part  of  the  Field  Data 
Collection  Program  which  is  managed  by  the  U.  S.  Army  Coastal  Engineer¬ 
ing  Research  Center.  The  U.  S.  Army  Engineer  Division,  South  Atlantic, 
and  the  U.  S.  Army  Engineer  Division,  New  England,  also  authorized  funds 
during  the  initial  year  of  this  study  (FY  1978)  to  expedite  execution 
of  the  Atlantic  coast  portion  of  this  program. 

This  report,  the  eleventh  in  a  series,  is  a  technical  report  de¬ 
scribing  an  integration  routine  to  calculate  nonlinear  energy  transfer 
for  waves  in  deep  water.  This  work  was  done  in  the  Hydraulics  Laboratory 
under  the  direction  of  Mr.  H.  B.  Simmons,  Chief  of  the  Hydraulics  Labora¬ 
tory,  Dr.  R.  W.  Whalin,  Chief  of  the  Wave  Dynamics  Division,  and 

Mr.  C.  E.  Chatham,  Jr.,  Chief  of  the  Wave  Processes  Branch.  This  report 

was  prepared  by  Mrs.  B.  A.  Tracy  and  Dr.  D.  T.  Resio.  Mrs.  D.  S.  Rags¬ 

dale  provided  computer  assistance,  and  Miss  C.  Lanford  prepared  some  of 
the  graphs. 

Commanders  and  Directors  of  WES  during  the  conduct  of  the  study 
and  the  preparation  and  publication  of  this  report  were  COL  John  L. 
Cannon,  CE,  COL  Nelson  P.  Conover,  CE,  and  COL  Tilford  C.  Creel,  CE. 
Technical  Director  was  Mr.  F.  R.  Brown. 
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THEORY  AND  CALCULATION  OF  THE  NONLINEAR  ENERGY 
TRANSFER  BETWEEN  SEA  WAVES  IN  DEEP  WATER 


Introduction 


1.  The  sea-wave  spectrum  usually  has  its  main  peak  just  before  a 
low  wave-number  cutoff.  The  JONSWAP  (Sell  and  Hasselmann  1972)  experi¬ 
ments  have  shown  that  this  peak  is  higher  and  narrower  than  had  been  ex¬ 
pected,  and  the  peak  moves  to  a  lower  wave-number  cutoff  as  time  in¬ 
creases.  Nonlinear  interactions  between  sea  waves  contribute  an  energy 
transfer  that  could  be  an  explanation  for  the  peak  enhancement  and  the 
shift  to  a  lower  wave-number  cutoff.  Integral  values  for  the  nonlinear 
energy  transfer  for  a  given  spectrum  can  be  used  to  show  how  this  spec¬ 
tral  situation  changes  over  a  period  of  time. 

2.  Nonlinear  wave-wave  interactions  become  especially  important 
in  the  case  of  tropical  wind  systems,  and  the  parameterization  of  wave- 
wave  interaction  among  components  of  an  asymmetric  spectrum  is  important 
for  the  case  of  tropical  systems  with  rapidly  varying  wind  fields.  Ex¬ 
perimental  data  (Forristall  et  al.  1978)  show  that  the  spreading  func¬ 
tion  is  frequency-dependent;  therefore,  different  frequency  components 
can  propagate  in  significantly  different  directions.  Investigators  feel 
that  the  nonlinear  transfer  may  be  the  reason  for  extremely  long  low- 
energy  waves  observed  after  severe  storms  (Hasselmann  1962).  An  integra¬ 
tion  technique  that  is  valid  for  the  nonlinear  interactions  in  a  sym¬ 
metric  spectr'im  can  be  easily  extended  to  an  asymmetric  spectral  system 
if  that  system  has  an  analytical  representation. 

3.  Hasselmann  (1962)  wrote  the  equations  and  evaluated  the 
boundary  conditions  that  model  the  behavior  of  a  system  of  waves  under¬ 
going  nonlinear  interactions.  These  equations  involve  the  rate  of 
change  of  mean  energy  (or  mean  action  density)  at  each  wave  number. 

(This  energy  flux  was  evaluated  by  using  a  perturbation  method.)  Hassel¬ 
mann  formulated  these  equations  by  considering  the  waves  as  particles 
and  by  considering  the  whole  process  as  particle  scattering.  Webb  (1978) 
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has  developed  an  integration  approach  that  evaluates  this  rate  of  change 
of  the  action  density  by  evaluating  the  resulting  Boltzmann  integrals. 
Development  of  a  special  type  of  grid  for  this  integration  scheme  allows 
a  fast  calculation  of  the  nonlinear  energy  transfer  for  a  whole  spectrum. 
This  grid  modification  consists  of  a  polar  grid  with  equal  angle  rays 
and  with  radial  rings  spaced  by  a  geometric  progression.  A  comparison 
of  the  results  of  this  method  to  Webb's  results  for  a  Pierson-Moskowitz 
spectrum  is  shown.  Results  for  various  other  symmetric  spectral 
situations  are  also  shown. 

Theoretical  Development 

4.  Webb  (1978)  has  taken  Hasselmann's  equations  and  written  the 
mean  action  density  equations  in  terms  of  a  transfer  function  which 
gives  the  rate  one  wave  is  scattered  into  another.  The  action  density 
form  of  Hasselmann's  equation  is* 

dk2  dk3  C(£r  *2,  £3,  k4)  *  6(k1  +  k2  -  i^3  - 
•  6(u)3  +  oj2  -  oj3  -  w4)  •  [n1n3(n4  -  n2>  +  n2n4^n3  -  n^]  (l) 

-V 

In  this  equation,  n^  is  the  action  density  at  wave  number  k^  ,  and 
is  the  angular  velocity  at  ki  .  6(...)  is  the  Dirac  delta  func¬ 

tion  and  C(...)  is  the  coupling  coefficient.  Energy  conservation  is 
taken  into  account  by  the  angular  velocity  delta  function,  and  momentum 
is  conserved  by  the  wave-number  delta  function.  The  process  involves 
a  set  of  four  waves  with  wave  numbers:  k^,  k2,  k3,  and  k^  .  Webb 
introduces  a  transfer  function  T(k3,k3)  where 


*  For  convenience,  symbols  and  unusual  abbreviations  are  listed  and 
defined  in  the  Notation  (Apppendix  B) . 
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dnl  f  ->•  -» 

=  Jdlc 3  T(k1>^3)  and 

TC^.fy  =  2  jjd^2  6 (1<1  +  -  k3  -  £  ) 


6(Ul  +  u2  -  u3  -  u4)  •  ©(l^  -  £j  -  |^  -  ^3|) 


[nin3(n4  -  n2)  +  n2n4(n3  -  n^]  (2) 


6(x)  =  1  if  x  >  0 
e(x)  =  0  if  x  <  0 


x  =  lk!  -  k4l  ~  lk!  -  k3i 


The  0(x)  function  determines  a  section  of  the  integral  which  is  not  de¬ 
fined  due  to  the  assumption  that  is  closer  to  lc3  than  k2  . 

5.  The  transfer  function  allows  us  to  think  of  the  process  as 
particle  scattering  where  we  consider  the  four  waves  as  four  particles 
with  their  momentum  related  to  their  wave  number.  The  transfer  function 
gives  the  rate  at  which  k3  is  scattered  into  wave  k^  .  The  density 
products  can  be  thought  of  as  a  diffusing  and  a  pumping  density  term 
where  n^n3(n4  -  n2)  is  the  pumping  term  and  n2n4(n3  -  n^)  is  the 
diffusive  term.  The  transfer  integral  can  be  thought  of  as  a  sum  of  a 
diffusing  transfer  integral  and  a  pumping  transfer  integral. 

6.  The  first  consideration  in  the  evaluation  of  the  diffusing 
and  pumping  transfer  integrals  is  the  limiting  properties  of  the  energy 
and  momentum  delta  functions.  We  will  use  the  following  property  of 
the  Dirac  delta  function  (Jackson  1962): 


6(x  -  a)  =  0  if  x  /  a 


where  x  and  a  are  representative  functions.  J 5(x  -  a)dx  =1  if 
region  includes  x  =  a  ,  and  is  zero  otherwise.  Consider 
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the  dk.  part  of  the  transfer  integral: 
4 


Jdit4  <5(1^  +  k2  -  k3  -  it4) 


If  +  k2  -  =  k4  for  conservation  of  momentum,  the  dk4  integral 

equals  one  and  the  transfer  integral  becomes 


T  (k^,  k^)  =  2jdk?6(u1  +  -  c.^)  Ffk^  ...) 


(4) 


where 


F  (k^ ,  ...)  =  CCi^,  .  .  .)e(x)  [n1n3(n4  -  n^  +  n2n4(n3  -  n.^] 

This  delta  function  evaluation  limits  the  wave  number  configuration  to 

— ►  — ►  — ^  «-V  — ►  — ► 

k,  +  k-  -  k_  =  k,  or  k,  -  k_  =  k,  -  k_  .  This  means  that  the  tips  of 

1/34134/ 

the  four  wave-number  vectors  must  form  a  parallelogram  in  wave-number 
space.  See  Figure  1  for  a  description  of  what  a  wave-number  configura¬ 
tion  would  look  like. 

7.  In  order  to  evaluate  the  integral  numerically,  we  must  fix 
values  for  k^(x,y)  and  l<3(x,y)  and  consider  the  limiting  properties 
of  the  angular  velocity  delta  function.  Then,  for  each  set  of  ( lc^ » k 3 )  , 
let  W(k2)  equal  the  argument  of  the  angular  velocity  delta  function. 

To  eliminate  the  angular  velocity  delta  function,  let  W(k2)  be  equal 
to  zero;  then 

W(k2)  =  0  =  +  w(£2)  -  o>3  -  u>4  (k^  +  k2  -  k3>  (5) 

The  momentum  conservation  condition  has  been  used  to  rewrite  k4  . 

-V  ^  ->  I* 

Consider  a  k^,  k^,  k3  coordinate  system — one  point  in  the  k^  -  k3 
plane  would  have  a  whole  line  of  solutions  parallel  to  the  k7-axis.  The 
set  of  solutions  that  will  satisfy  the  conservation  conditions  can  be 
represented  as  an  egg-shaped  two-dimensional  locus  in  a  Cartesian  coor¬ 
dinate  system  in  I<2  space  where  k^  is  the  x-axis  and  k^  is  the 

y-axis.  On  this  locus  n  is  in  the  normal  or  radial  direction  and  s 
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Figure  2.  A  (n,s)  coordinate  system  is  set  up  using 
the  normal  and  tangential  vectors  relative  to  the  locus 


Another  property  of  the  Dirac  delta  function  is  (Jackson  1962): 


6  (x  - 


6 [ f (x) 


x  ) 
o 


djf 

dx 


where  f(x  ) 
o 


[  6 ( f ) d f  =  6(x)dx ] 


where  x  and  x  are  representative  functions.  Using  this  property  on 

->  -v  ->  ->•  (5  ( —  0) 

dn • [W(s,n) ]  ,  we  find  that  6[W(s,n)]  =  — - —  .  We  must  integrate 


around  the  s-curve  for  this  property  to  hold,  so 


TO^.tj)  =2  ff ds  dn  •  F(k1,  ...) 

9W(s,n) 

-  ^ f-~-  f<ki.  •••>  <7> 

9W(s,n) 

3  a 


since  Jdn  S(n  -  0)  =1  .  (Note  that  evaluation  of  the  integral  now  re¬ 
duces  to  evaluation  of  a  line  integral  or  contour.)  The  normal  deriva¬ 
tive  in  the  denominator  is  the  magnitude  of  a  gradient  and  can  be 
written : 


3VJ(s  ,n) 
3n 


-1 


=  l™<k2x,k2y) 


i-l 


(8) 


This  term  can  be  thought  of  as  a  phase  space  term. 

9.  The  evaluation  of  the  integral  now  entails  evaluation  of  all 
the  factors  of  the  integrand  at  a  series  of  evenly  spaced  points  along 
the  k^-contour  loop  or  locus.  The  factors  of  the  integrand  include  the 
phase  space  term,  the  coupling  coefficient,  and  the  density  function. 
Each  of  these  will  be  dealt  with  separately.  Figures  3  and  4  show  how 
these  three  terms  behave  over  a  representative  locus.  A  representative 
locus  is  also  shown  in  Figure  5.  Fifty  evenly  spaced  increments  are 
shown  as  the  points  on  the  locus. 


10. 


The  phase  space  term 


DW 

'n 


can  be  evaluated  at  each  point  of 


the  locus  by  evaluating  the  gradient  of  the  locus  equation  (Equation  5). 

_  .  „  ,1/2 

For  deep  water  w  a  k 


1/2 


fied  by  using  Q  =  k1  -  k 


(k  =  (k|).  The  locus  equation  can  be  simpli- 

1/2 


(since  k^  and  k.^  will  remain  con¬ 


stant  during  calculation  of  the  k^-contour)  and  defining  a  new  P-vector, 

— i ►  — V  k 

P  =  k^  -  k^  .  Each  interaction  can  be  defined  by  a  specific  P-vector. 

~y 

After  evaluating  the  w- funct ions  and  substituting  for  Q  and  P  ,  the 


expression  for  Wfk^)  becomes: 


Q  +  k 


1/2 


(r?  +  s2i) 


1/2 


=  0 


(9) 
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Figure  3.  Behavior  of  the  phase  space  term  and  density  functions 
over  a  representative  locus 


Figure  4.  Behavior  of  the  coupling  coefficient  over  a 
representative  locus 


The  gradient  of  W  is: 


where 


and 


VW 


3W 

9x 


3H  a[k21/2-  (l?  +  t2|)1/2] 

8x  3x 


w  %n  -  (if  --  F2l)1/2] 

3y  3y 


(10) 


(ID 


(12) 
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LOCUS  FOR  k2 


Figure  5.  A  representative  locus  in  Cartesian 
coordinates.  Locus  coordinates  have  been 
found  by  an  iteration  routine  utilizing  binary 

search 


Writing  the  vectors  in  rectangular  components  yields: 


h +  *>2 


7-|lM 

+  (py  +  y)  J 


(13) 


where  P  and  P  are  the  rectangular  components  of  P  and  x  and  y 

x  y 

are  the  rectangular  components  of  .  The  magnitude  of  the  gradient 

in  the  normal  direction  is: 
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vw 


[(f)2 + (f  )2]1/2 


2  “3/4  (Px  +  x)  r  2  2 

,f  (X  +  y  )  -  — -  [(Px  +  x)2  +  <Py  +  y)2 


rf 


+  (*2  +  y2) 


-3/4  (P  +  y) 


[<p*  +  ;': 


.  1/2 


1-3/4# 


+  (P  +  y)‘ 


(14) 


For  ease  in  working  with  the  angles  involved,  rectangular  coordinates 
were  used  for  evaluation  of  this  term. 

11.  The  coupling  coefficient  from  Hasselmann's  work  is  a  complex 
algebraic  expression  and  in  Webb's  paper  was  put  in  a  usable  form  by  an 
algebraic  manipulator.  The  expression  given  in  Webb's  work  was  used, 
and  the  behavior  of  this  term  when  plotted  over  a  representative  locus 
exhibits  the  same  shape  as  Webb's  curve.  The  complete  algebraic  expres¬ 
sion  is  listed  in  the  appendix. 

12.  The  density  term  in  Hasselmann's  equation,  n,n_(n,  -  n_)  + 

i  3  4  2 

n^n^(n^  -  n^)  ,  can  be  thought  of  as  the  sum  of  two  types  of  transfer. 

A  diffusive  transfer  between  k^  and  k^  is  described  by  n^n^Cn^  -  n^) 

This  diffusive  transfer  causes  a  pumped  transfer  between  k^  and  k^ 

in  n,n_.(n.  -  n_)  .  The  n.  values  are  action  densities  and  are  re- 
13  4  2  i 

lated  to  the  amount  of  energy  inherent  in  a  given  state,  k^  .  The  action 
density  subscripts  correspond  to  the  wave-number  subscripts.  The  whole 
purpose  in  constructing  an  evaluation  scheme  for  the  Boltzmann-type 
integrals  is  to  evaluate  the  rate  of  change  of  the  action  density  at  a 
certain  wave  number  and,  therefore,  get  a  value  for  the  energy  transfer 
there.  Right  now,  we  are  interested  in  knowing  the  amount  of  action 
density  we  have  available  to  cause  this  diffusing  and  pumping  effect  in 
the  set  of  four  wave  numbers,  k^,  k^ ,  k^,  and  k^  . 

13.  The  n^  (k^)  values  are  directly  calculated  from  the  spectrum 


13 


that  is  being  consider?'4.  Most  spectra  can  be  written  as  an  E(f)  func- 

2 

tion  where  E  is  energy  in  in  /Hz  .  The  independent  variable  is  fre¬ 
quency,  but  we  are  dealing  in  wave-number  space  so  we  must  transform 
the  energy  density  function  from  a  function  of  frequency  to  a  function 

of  wave  number.  The  one-dimensional  frequency  function  was  transformed 

2 

to  a  two-dimensional  function  by  using  a  cos  fi  spreading  function. 

The  group  velocity  c  can  be  written  as: 


c 

g 


(15) 


where  w  is  angular  velocity.  The  identity  relating  frequency  to  wave 
number  in  deep  water  is 


u>  a  k  or  w  =  /gk 


(16) 


Then, 


_ w  /gk  _  _1 


2k 


2k 


(17) 


The  density  function  n^k^  can  be  written  as: 


n1(ki)  = 


F(k.) 


(18) 


where  F(k^)  is  the  two-dimensional  spectrum  with  respect  to  wave  num¬ 


ber  and  o)^  is  the  angular  frequency.  The  c^ 


term  acts  like  a  modula¬ 
tion  term  on  the  amplitude  of  the  wave  function  making  up  the  E(f) 
term  to  enable  a  transformation  from  a  frequency-dependent  function  to  a 
wave-number  dependent  function.  The  whole  procedure  amounts  to  multi¬ 
plying  E(f)  by  the  square  of  c  /w  and  adding  the  spreading  function 

g  2  2 

by  multiplying  the  whole  expression  by  —  •  cos  0  .  If  we  know  the  wave 
number,  the  action  density  at  a  specific  wave  number  can  be  calculated 
using  the  above  relationships.  A  Pierson-Moskowitz  spectral  density  func¬ 
tion  was  used  in  the  initial  test  case  to  match  with  Webb's  work,  but  any 
spectrum  with  an  analytical  wave-number  representation  could  be  used  in 
the  calculations. 
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14.  So  far,  the  basic  aspects  of  the  factors  of  the  integrand  have 
been  dealt  with.  All  these  hinge  on  calculation  of  the  k„-contour  or 
k^-locus  for  each  k^  -  k^  interaction.  It  is  possible  to  calculate  a 
locus  by  iteration  or  a  combined  interpolation-iteration  scheme  for 

each  separate  (k^,k^)  couplet,  and  this  would  be  the  approach  to  the 
problem  if  we  used  a  regularly  spaced  (k,  0)  grid.  However,  using  a 
geometric  progression  spacing  scheme  for  the  k  (radial)  values  allows 
us  to  utilize  some  geometric  vector  relationships  that  will  permit  the 
estimation  of  the  k2~  and  k^-loci  without  having  to  calculate  each 

-y  ~y 

as  an  independent  function  of  a  (k^,k^)  couplet  after  one  locus  has 

been  iterated  for  each  possible  orientation  of  the  (k^,k^)  couplet. 

In  this  method  an  initial  value  is  chosen  for  the  magnitude  of  k^ 

where  k^  is  the  smallest  value  used  for  r  on  the  (r,  0)  grid.  We 

defined  our  value  for  k  to  be  0.14  m  ^  (MKS  wave-number  units)  . 

o 

k  must  be  chosen  to  be  a  small  wave  number  in  the  area  where  very 
o 

little  interaction  takes  place.  In  k-space  we  are  using  polar  coordi¬ 
nates  where  |k|  =  r  and  0  is  incremented  from  0°  to  90°  in  9°  incre¬ 
ments.  Our  first  area  of  reference  will  be  a  circle  of  radius  =  k 

o 

The  next  reference  circle  around  the  origin  will  have  r  =  Xk  ,  where 

°  2 

we  have  specified  X  =  1.2  .  The  next  concentric  circle  will  be  X  kQ 

and  so  on.  Figure  6  shows  the  completed  grid  with  X  =  1.2  and 

k  =  0.14  . 
o 

15.  For  each  interaction  we  need  to  define  a  ^-vector 

(P  =  k^  -  lc^)  •  This  P-vector  will  define  the  orientation  of  the  inter¬ 
action  on  the  grid.  As  stated  before,  each  interaction  consists  of  a 
constant  k^-vector  and  a  constant  k^-vector.  | k^ |  will  be  assumed  to 
be  greater  than  |ic^|  for  ease  in  calculation.  ?  can  be  calculated 
from  the  k^-  and  k^-vectors,  and  there  exists  a  specific  ?-vector  for 
each  interaction.  Considering  Figure  7,  it  is  obvious  that  on  rays  of 
constant  0  a  set  of  parallel  ^-vectors  is  produced.  These  parallel 
^-vectors  can  be  related  to  the  locus  of  the  initial  ^-vector  by  the  X 
scaling  factor. 


ORIGIN 


iW!\ : 


tiV’AV 


Figure  7.  The  geometric  construction  of  a  set  of  parallel 

interaction  vectors 


16.  Considering  the  geometry  in  Figure  7  with  sin  6^  =  y^/|xk£ 

for  triangle  OBC,  the  Cartesian  coordinates  for  the  orientation  of 
“► 

P.  (one  of  the  interaction  P-vectors  that  utilize  the  k  -vector)  are: 
J-  o 


y  =  X  k  sin  6. 

o  1 


x  =  A  k  cos  6.  -  k 
o  1  o 


then. 


The  components  of  can  be  evaluated  in  the  same  way  using  triangle 

ODE  in  Figure  7 : 


Looking  at  the  ratio  of  j  |  to  |p  |  ,  we  find  | |  /|p  |  =  X  .  So, 
the  ratio  of  the  magnitude  of  two  consecutive  parallel  P-veclors  is  X  . 

17.  Considering  the  locus  equation  and  the  relationship  between 
|?2 I  and  |?^|  ,  we  can  show  that  if  we  have  the  locus  for  the  P^ 
interaction,  we  can  calculate  the  locus  for  the  P2  interaction.  The 
locus  equation  is: 

Q  +  k*/2  ~  (|P  +  k2l)1/2  =  0  (21) 

where 
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Q  =  k2/2  -  k2/2  and  ?  =  -  ic, 


First,  let  us  write  the  k^  and  k^  interaction  components  (Figure  8) 
for  the  locus  equation  for  : 


k.  =  Xk  :  k  =  Xk  cos  9,  k  =  Xk  sin 
3  o’  x^  o  y3  ° 


k,  =  k  ;  k  =k,k  =0 
1  0  x!  °  yj_ 


(22) 

(23) 


Second,  the  k^  and  k^  interaction  components  for  the  locus  for 
|?2|  =  X  I, P 1 1  are: 


Xk  ;  k  =  Xk  ,  k  =  0 
ox,  o  y. 


k„  =  X2k  ;  k  =  X2k  cos  9,  k  =  X2k  sin  6 
3  0  x3  °  y3  o 


(24) 

(25) 


Substituting  the  above  information  we  find  for  the  first  locus  equation: 

Q  =  k^/2  -  (XkQ)1/2  =  k2/2(l  -  X1/2) 


P,  =  k  -  Xk  cos  0 
lx  o  o 


P,  =  0  -  Xk  sin  0 
ly  o 


1*1  +  4l  =  >Plx^k2x)2  +  (Ply  +  V' 


By  algebraic  manipulations  and  substitution. 


I*5!  +  "  yjk2  +  k2x(2ko)(1  "  ■*  cos  0)  +  k2y(2>lt0>(-sin  8)  +  k^(l  +  \2  -  2\  cos  B)  (^6) 


Using  the  general  locus  equation, 

q  ♦  k!/2  - 


(lf^2l) 


1/2 


=  0 


(27) 
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the  first  locus  equation  becomes 


k*/2(l  -  X1/2)  +  k^/2  -  £k2  +  k2x(2ko)(l  -  X  cos  0) 


2  2  r/4 

+  k  (2Xk  ) (-sin  0)  +  k  (1  +  X  -  2X  cos  6)  =  0 

2v  o  o  J 

Substituting  the  information  for  the  second  locus,  we  find 

q  -  Ok  >1/2  -  (A  )1/2  .  k1/2J1/2o  -  >1/2> 

o  \  0/  o 


(28) 


Pn  =  Xk  -  X  k  cos  0 
2xo  o 


P0  =  0  -  Xk  sin  0 
2y  o 


l?2  +S2'  ’  V(P 2*  +  k2,>2  + 


(P2y  +  k2y} 


(29) 


By  algebraic  manipulation  and  substitution. 


l?2  +  k2l  k2  +  k2x(2koX)(1  "  '  cos  9)  +  k2  <-2ko*'2  S*n  6)  +  '2ko(>2  '  2'  l'OS  6  +  (30) 


The  second  locus  equation  becomes 


k1/2d/2 

o 


,,  .1/2.  ,  .1/2 
(1  -  X  )  +  k2 

k2  +  k„  (2k  A)(l  -  X  cos  6) 

1  l  Zx  0 

+  k  (-2k  X2  sin 
2y  V  o 

\  2  2  2  "Jl/^ 

9j  +  A“V(A  -  2X  cos  0  +  1)  =0 

(31) 

18.  Using  the  assumption  that  the  loci  scale  by  a  factor  of  X, 


we  should  be  able  to  substitute  k„  =  Xk„  ,  k„  =  Xk„  ,  and 

n - 2 - 2 - 2  2X2  2x1  2y2  2yi 

k.  =  .  /X  k„  +  X  k„  into  the  second  locus  and  have  an  algebraic  ex- 
2  \  2x^  2yx 

pression  equivalent  to  the  expression  for  the  first  locus.  Substituting, 


the  second  locus  becomes. 


1/21/2  1/?  /  2  2  9  9 

k  A  7  (1  -  A  '~)  +  1  *  Z 

o 


(2  2  9  2 


2  2  2  2 
X  k,  +  >.  k, 

2xl  2yl 


+  k„  (2k  X)(l  -  X  cos  9)  +  k.  (-2k  X^  sin  0) 
2xo  2v  o 


+  >2k2(X2  -  2)  cos  0+1) 
o 


1/4 


0  (32) 


1/2 

Divide  the  above  equation  by  X  and  it  is  equivalent  to  the  expres¬ 

sion  for  the  first  locus  equation.  Therefore,  if  we  have  the  coordinates 

of  the  initial  locus,  we  can  calculate  the  locus  coordinates  for  any 

— ►  — ► 

interaction  that  has  a  P-vector  parallel  to  the  P-vector  of  the  initial 
interaction. 

19.  In  calculation  of  the  loci,  we  need  only  calculate  a  basic  set 

of  loci  where  k^  =  ,  9^  =  0°  ,  and  k^  takes  the  position  of  each 

other  point  of  intersection  of  the  grid.  Comparisons  are  shown  in 

Table  1.  In  this  figure  we  have  considered  a  basic  locus  where 

k,  =  0.1400  ,  k,  =  0  ,  k„  =  0.2016  ,  and  k„  =  0  .  The  x  and  v 
lx  ly  3x  3y 

coordinates  are  given  for  this  locus.  The  next  set  of  values  are  X 
times  the  coordinates  of  the  basic  locus.  The  third  section  contains 
the  coordinates  for  a  locus  with  the  same  P-vector  coordinates  as  the 
second  locus.  The  basic  locus  and  the  locus  of  the  third  section  have 
been  iterated  by  a  binary  search  routine  that  will  be  outlined  later  in 
the  section  on  the  program  method.  (k^  =  0.1680  ,  k^  =  0  , 
k^x  =  0.2419  ,  and  k^  =  0  are  the  interaction  components  of  the 
second  and  third  loci.)  The  ratio  of  the  P-vector  of  the  second  locus 
to  the  ^-vector  of  the  basic  locus  is  X  .  It  is  obvious  that  X  times 
the  basic  locus  coordinates  gives  locus  coordinates  that  are  very  close 
to  the  iterated  locus  values.  Figures  9  and  10  and  Table  2  show  two 
loci  where  the  ratio  of  P^  to  is  1.15  . 

20.  The  solution  to  the  problem  of  the  calculation  of  the  loci 

reduces  to  the  calculation  of  a  basic  set  of  loci.  The  loci  for 
— ► 

parallel  P-vectors  'or  other  higher  order  interactions  can  be  calculated 
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Table  1 

Comparison  ol  a  Basic  Locus  and  a  Scaled  Locus 


Basic  Locus 

tor 

k,  =  0.1400  ,  k,  =  0.0000 
lx  ly 

k  =  0.2016  ,  k,  =  0.0000 

lx  ly 

*  “  1  .  i. 

Lambda  Times 

Basic  Locus  Values 

Iterated 

for 

k.  =  0.1680  , 
lx 

k  =  0.2419  , 
3x 

Locus 

k.  =  0.0000 

ly 

k  =  0.0000 

3y 

X 

y 

X 

y 

X 

V 

0.101b 

0.0 

0.2419 

0.0 

0.2419 

0.0 

0.2007 

0.0098 

0.2408 

0.0118 

0.2409 

0.0118 

0.1982 

0.0193 

0.2378 

0.0232 

0.2378 

0.0232 

0.1920 

0.0282 

0.2328 

0.0335 

0.2329 

0.0338 

0.1886 

0.0362 

0.2263 

0.0433 

0.2264 

0.04  34 

0.1822 

0.0432 

0.2186 

0.0517 

0.2187 

0.0517 

0.1750 

0.0489 

0.2100 

0.0587 

0.2100 

0.0587 

0.1674 

0.05  37 

0.2009 

0.0644 

0.2008 

0.0645 

0.1594 

0.0575 

0.1913 

0.0690 

0.1913 

0.0690 

0.1511 

0.0603 

0.1816 

0.0724 

0.1815 

0.0723 

0.1431 

0.0622 

0.1717 

0.0746 

0.1717 

0.0746 

0.1350 

0.0633 

0.1620 

0.0760 

0.1620 

0.0759 

0.1269 

0. 0636 

0.1523 

0.0763 

0.1523 

0.0763 

0.1189 

0.0633 

0.1427 

0.0760 

0.1427 

0.0759 

0.1110 

0.0625 

0.1332 

0.0748 

0.1332 

0.0747 

0 . 1032 

0.0606 

0.1238 

0.0727 

0.1239 

0.0727 

0.0955 

0.0583 

0.1146 

0.0700 

0.1146 

0.0699 

0.0879 

0.0523 

0.1055 

0.0662 

0.1054 

0.0663 

0.G803 

0.0515 

0.0964 

0.0618 

0.0964 

0.0618 

0.0730 

0.0469 

0.0876 

0.0563 

0.0876 

0.0563 

0.0659 

0.0415 

0.0791 

0.0498 

0.0790 

0.0497 

0.0592 

0.0340 

0.0710 

0.0420 

0.0710 

0.0421 

0.0532 

0.0276 

0.0638 

0.0331 

0.0638 

0.0331 

0.0482 

0.0192 

0.0578 

0.0230 

0.0579 

0.0230 

0.0448 

0.0099 

0.0538 

0.0119 

0.0538 

0.0118 

0.0436 

0.0000 

0.0523 

0.0000 

0.0524 

0.0000 

0.0448 

-0.0099 

0.0538 

-0.0119 

0.0538 

-0.0118 

0.0482 

-0.0192 

0.0578 

-0.0230 

0.0579 

-0.0230 

0.0532 

-0.0276 

0.0638 

-0.0331 

0.0635 

-0.0331 

0.0592 

-0.0350 

0.0710 

-0.0420 

0.0710 

-0.0421 

0.0659 

-0.0415 

0.0791 

-0.0498 

0.7090 

-0.0497 

0.0730 

-0.0469 

0.0876 

-0.0562 

0.0876 

-0.0563 

0.0803 

-0.0513 

0.0964 

-0.0618 

0.0964 

-0.0618 

0.0879 

-0.0552 

0.1055 

-0.0662 

0.1054 

-0.0663 

0.0955 

-0.0582 

0.1146 

-0.0698 

0.1146 

-0.0699 

0.1032 

-0.0606 

0.1238 

-0.0727 

0.1239 

-0.0727 

0.1110 

-0.0623 

0.1332 

-0.0748 

0.1332 

-0.0747 

0.1189 

-0.0633 

0.1427 

-0.0760 

0.1427 

-0.0759 

0.1269 

-0.0636 

0.1523 

-0.0763 

0.1523 

-0.0763 

0.1350 

-0.0633 

0.1620 

-0.0760 

0.1620 

-0.0759 

0.1431 

-0.0622 

0.1717 

-0.0746 

0.1717 

-0.0746 

0.1513 

-0.0603 

0.1816 

-0.0724 

0.1815 

-0.0723 

0.1594 

-0.0575 

0.1913 

-0.0690 

0.1913 

-0.0690 

0.1674 

-0.0537 

0.2009 

-0.0644 

0.2008 

-0.0645 

0.1750 

-0.0489 

0.2100 

-0.0587 

0.2100 

-0.0587 

0.1822 

-0.0431 

0.2186 

-0.0517 

0.2187 

-0.0517 

0.1886 

-0.0361 

0.2263 

-0.0433 

0.2264 

-0.0434 

0.1940 

-0.0282 

0.2325 

-0.0338 

0.2329 

-0.0338 

0.1982 

-0.0193 

0.2378 

-0.0232 

0.2378 

-0.0232 

0.2007 

-0.0098 

0.2408 

-0.0118 

0.2409 

-0.0118 

Table  2 


A  Comparison  of  the  Two  Loci  in  Figures  9  and  10 


LOCUS  1 

LOCUS  2 

RATIOS 

; 

k,  =0.100 

k,  =0.115 

UM>2 
<Mi  ' 

1.1500 

k3=0.115 

k3  =  0.133 

(k3)2  _ 

1 

•  1.1500 

P,  =0.0575 

P2  =  0.066 

<P>2 
<P)l  = 

1.1478 

X 

y 

X 

y 

■1 

0.5756 

0.7806 

0.6636 

0.9000 

1.1C29 

1.1530 

0.3891 

1.3550 

0.4486 

1.5623 

1.1529 

1.1530 

-0.1155 

1.5151 

-0.1334 

1.7477 

1.1550 

1.1535 

-0.4703 

1.2453 

-0.5425 

1.4360 

1.1535 

1.1535 

-0.5763 

0.8565 

-0.6648 

0.9875 

1.1536 

1.1529 

-0.5032 

0.4904 

-0.5805 

0.5653 

1.1536 

1.1527 

-0.2664 

0.1622 

-0.3072 

0.1871 

1.1532 

1.1535 

0.1573 

0.0852 

0.1813 

0.0982 

1.1526 

1.1526 

0.4494 

0.3817 

0.5188 

0.4394 

1.1544 

1.1512 

0.5676 

0.7120 

0.6547 

0.8208 

1.1535 

1.1528 
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from  the  basic  set.  Two  other  factors  of  the  integrand. 


aw 

3n 


and  the 


coupling  coefficient,  also  scale  by  functions  of  the  A-scaling  factor 
for  the  parallel  P-vectors. 


DW  =  |3W/3n 


21.  The  phase  space  term, 
-1 


-1 


,  scales  by  the  /\  .  Let 


for  the  basic  locus  for  the  interaction  listed  in 
-1 


Table  1.  Let  DW^  =  |3W/3n|  *  for  the  locus  that  has  a  P-vector 
parallel  to  the  P-vector  which  is  A  times  the  first  P-vector  as 
shown  in  the  second  column  of  Table  1.  We  find  that  DW2/DW^  =  /a  . 
These  ratios  are  shown  in  Table  3.  The  eight  values  listed  here  are 


the  calculated 


3W/3n 


i-l 


values  for  the  first  eight  locus  coordinates 


shown  in  Table  1.  Similarly,  the  coupling  coefficient  scales  by  a 

factor  of  A*3  .  Values  of  the  calculated  coupling  coefficient  for  two 

loci  are  shown  in  Table  4  for  the  same  locus  coordinates  as  were  used 

in  the  above  discussion  on  the  phase  space  term. 

22.  The  | 3W/ 3n |  ^  term  as  stated  before  is  the  magnitude  of  the 

/  2  2  \  1  /  2 

gradient  of  the  W  expression  for  the  locus  where  =  fk^  +  k^  1 


VW 


■  M  *  (f /] 


1/2 


'  (|i  (“*  +  O' 


k  -3/4  (P  +  k  )  r 

U  +  kz)  — *■-?-*  up  +  k  r 

x  y)  2  l  x  x 


+  (P  +  k  ) 

y  y 


■r  ’  •  14  « *  .r" 


(33) 


+  ^y^  I  2  2 

- *—  l(P  +  k  )  +  (P  +  k 

xx  y  y 


[ 


-3/4 


Vl/2 


Let  k  =  A  k 
x  x 


and  k  =  A  k  and  P  =  AP  .  Then 

y  y 
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Table  3 

Scaling  of  the  Phase  Space  Term 


>  Phase  Space  Term 
,jW;  -1 
:  bl 

t 

Phase  Space  Term 

I  0  W'  “1 

!«;  »BL 

Ratio  of 

Phase  Space  Terms 

,3W,  -1  ,,3Wi  -1 

l5n'  XBL  3n'  BL 

Basic  Locus 

a  Times  Basic  Locus 

X  =  1.2 

X  =  1.2 

o 

o 

<T 

H 

r— 

kx  =  (0.168, .00) 

VT  =  1.095 

=  (0.2016, .00) 

k,  =  (0.242, .00) 

3 

4.49 

4.92 

1.096 

4.44 

4.87 

1.097 

4.31 

4.  72 

1.095 

4.11 

4.  50 

1.095 

3.85 

4.22 

1.096 

3.58 

3.92 

1.095 

3.29 

3.60 

1.094 

3.00 

- - - - - - - - 

3.29 

1.097 

Table  4 

Coupling  Coefficient 


Coupling  Coefficient 

CBL  X  10'3 

Coupling  Coefficient 

CXBL  *  10"3 

Ratio  of 

Coupling  Coefficients 

CXBL/CBL 

Basic  Locus 

X  Times  Basic  Locus 

X  =  1.2 

X  -  1.2 

k  =  (0.14, .00) 

k3  =  (0.2016,. 00) 

kx  -  (0.168,. 00) 
k3  -  (0.242,. 00) 

X6  -  2.986 

0.196 

0.586 

2.990 

.193 

.576 

2.984 

.183 

.547 

2.989 

.168 

.502 

2.988 

.149 

.446 

2.993 

.129 

.384 

2.976 

.108 

.323 

2.991 

.0885 

.264 

2.983 
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(34) 


23.  All  these  scaling  factors  will  aid  in  the  numerical  evalua¬ 
tion  of  the  transfer  integral  by  allowing  us  to  use  one  basic  set  of 
calculations  and  the  scaling  factors  to  evaluate  the  interactions  for 
the  entire  grid.  The  transfer  integral  can  be  written  as: 

T  (k^  £3)  =  f  F  (k^  ...)  D  (kr  ...)  (35) 


28 


where 


F  (k1,  ...)«-§*  *  C(k1,  k2,  k3,  k4)  •  ds 

Density,  D  (k  ,  ...)  ,  is  the  sum  of  the  diffusing  and  pumping  densities 
for  a  given  set  of  four  wave  numbers  and  a  given  spectral  representation. 
Note  that  the  0(x)  function  has  been  left  out  of  the  integral  repre¬ 
sentation  since  we  will  not  delete  a  section  of  the  contour  as  Webb  has. 
The  transfer  integral  can  be  written  as  a  product  of  a  geometric  product 
term  and  a  density  function.  The  density  function  is 

nln3 (n4  "  n2^  +  n2n4^n3  "  ni^  (36) 

where  n.  is  the  action  density  at  various  wave  numbers,  k.  .  This 
term  can  be  calculated  using  the  locus  coordinates  for  the  specific 
interaction  and  the  applicable  spectral  function.  The  F  (k^ ,  ...) 
term  is  a  geometric  product  term  depending  on  the  specific  locus.  When 
the  basic  locus  coordinates  are  calculated,  a  geometric  product  term  is 
also  calculated  for  each  increment  of  the  basic  locus.  This  geometric 
product  factor  can  be  scaled  for  A-multiples  of  the  basic  loci  by  a 
product  of  the  separate  scaling  terms  as  defined  before.  This  product 
of  scaling  terms  is 

/l  •  A6  •  A  =  A15/2  (37) 

The  last  A  is  obtained  from  the  polar  differential  increment  for  the 
integral,  r  dr  d0  or  ds  .  The  transfer  integral  can  then  be  written 
as  T($1,^3>  =  Jd  (k^,  ...)  •  GEOM  (k^,  ...)  •  A15^2  where  the  A 
value  is  determined  by  the  interaction  being  considered.  The  basic  geo¬ 
metric  product,  GEOM  (k  ,  ...)  ,  is  determined  by  the  basic  locus  that 
has  a  ^-vector  orientation  parallel  to  the  P-vector  of  the  interaction. 

24.  The  integral  in  the  above  form  has  advantages  when  using  com¬ 
puter  evaluation.  The  product  form  has  a  vectorizing  potential  that 
will  make  evaluations  of  the  nonlinear  energy  transfer  for  a  whole  spec¬ 
trum  fast  and  inexpensive. 

25.  Going  back  to  the  theoretical  development,  the 
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nonlinear  energy  transfer  can  be  written  as: 


3iij 

W 


fu^L,i3)  dk. 


(38) 


In  polar  coordinates 

T(it1,k3>  de3k3dk3  (39) 

As  discussed  before,  the  transfer  at  a  specific  orientation  of  k^  and 
k3  will  have  an  initial  value  which  is  a  function  of  the  basic  locus. 

-t  4 

This  basic  value  of  the  transfer  integral  will  be  defined  as  To(k3,k3) 
Each  successive  transfer  integral  at  the  same  orientation  of  k^  and  k 
(parallel  P-vectors)  can  be  written  as  A^^  •  T 

o 

26.  Therefore,  at  a  fixed  k^  -  k3  orientation  (0^)  ’  t^ie  non~ 
linear  energy  transfer  can  be  written: 

3ni  7  n  1S/9 

“  (013  fixed)  =  /  (A0)1  '  To  k3  dk3  (40) 

Jo 

n+1 

k3  can  be  written  as  A  kQ  for  a  specific  interaction  where  n 
(the  number  of  the  radial)  begins  at  zero  and  progresses  to  infinity  (in 
reality,  the  edge  of  the  grid).  dk3  or  Ak3  can  be  written  as 
(An+^  -  An)  kQ  .  In  summation  notation,  the  transfer  for  one  angular 
section  is: 


3n 


~  (013  fixed)  =  X  (3n)15/2  Tq  An+1  ku(An+1  -  An)kQ  (41) 
n=0 


where  n  is  the  number  of  the  radial  and  infinity  is  the  edge  of  the 
grid.  The  total  nonlinear  transfer  for  the  whole  grid  can  be  written  as 


3nl  V  3ni 

3t  ^  at  ^ei3  fixed) 


A0 


13 


(42) 


where  the  summation  is  over  all  possible  0^3  orientations. 
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Program  Method 

27.  A  polar  geometrically  spaced  grid  as  described  in  the  pre¬ 

vious  section  (Figure  6)  was  set  up  and  used.  The  numerical  results 
given  in  this  report  were  obtained  by  using  a  grid  with  =  0.14  in 

MKS  wave-number  units  and  >.  =  1.1  .  The  6  increments  on  the  grid 
were  4.5°.  Thirty  rings  of  different  r  values  are  used  and  eleven 
positive  and  negative  angle  increments  including  zero  were  used. 

28.  A  set  of  basic  loci  was  calculated  using  the  above  grid 
parameters.  This  was  done  in  a  separate  program,  and  the  output  of  this 
program  created  a  tape  listing  the  basic  locus  coordinates  with  their 
corresponding  geometric  product  calculation  for  each  possible  interaction 
angle.  To  determine  the  basic  loci  k^  assumed  the  position  of  each 
intersection  on  the  grid  and  k^  remained  equal  to  k^  and  oriented 

at  0=0°.  All  the  angular  points  of  each  ring  were  done  first,  and 
then  the  program  moved  on  to  the  next  ring.  The  loci  were  stored  in 
this  order. 

29.  To  locate  the  coordinates  for  each  of  the  basic  loci  the  ana¬ 
lytical  expression  (Equation  9)  for  the  locus  equation  was  solved  using 
k^  and  k^  to  find  the  coordinates  (k^)  of  a  point  on  either  end  of 
the  egg-shaped  locus;  then  the  center  for  the  locus  was  found  by  locating 
the  middle  point  of  the  line  between  two  points  on  the  ends  of  the  locus. 
A  radius  for  the  locus  was  calculated  from  this  information  and  a  binary 
search  procedure  found  50  evenly  spaced  points  of  the  locus  within  an 
epsilon  range  of  10  The  first  coordinate  of  each  locus  is  the  point 
where  the  P-vector  of  the  interaction  would  intersect  it.  The  points 
continue  in  a  counterclockwise  direction  until  the  loop  is  complete. 

All  the  coefficients  of  the  integrand  (the  coup  ing  coefficient,  the 
phase  term,  and  the  ds  term)  except  the  densities  are  combined  into 
a  product  and  saved  for  each  locus  increment.  Each  of  these  integrand 
coefficients  was  calculated  via  a  subroutine  for  each  point  on  each  locus. 
These  basic  locus  coordinates  and  a  value  for  the  product  (geometric 
product  factor)  of  the  integrand  coefficients  are  available  to  be  read 
into  the  main  program.  The  basic  loci  were  read  into  the  main  program 
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and  indexed  by  the  difference  of  the  k^-angle  and  the  k^-angle,  the  ring 

difference  between  and  k^  ,  and  the  number  of  the  increment  on  the 

specific  locus  for  the  k1  and  interaction.  The  geometric  products 

were  also  indexed  by  this  same  scheme.  Loops  were  set  up  to  define  the 

| k |  values  for  each  point  of  the  grid  using  the  >-scaling  factor  and 

15/2 

to  define  the  value  of  the  geometric  scaling  factor  (,\  )  for  each 

point  of  the  grid.  Loops  were  also  set  up  to  determine  the  values  for 
the  angles  used  in  the  calculation. 

30.  The  incrementation  scheme  in  the  integration  program  uses  a 

series  of  four  nested  do  loops.  The  outermost  loop  is  the  angle  of 

k^  .  Directly  inside  this  loop  we  increment  k^  by  the  ring.  The 

location  of  k^  and  the  angle  of  k^  will  determine  the  position  of 
— ► 

the  P-vector  for  the  interaction.  The  ring  location  of  k^  will  deter¬ 
mine  the  A-scaling  factor  for  the  corresponding  set  of  k^  -  k^  inter¬ 
actions  that  have  the  same  angle  orientation  as  .  This  will  allow 

us  to  calculate  the  locus  coordinates  and  the  geometric  product  factors 
for  all  the  interactions  in  this  set.  The  ^-incrementation  takes  place 
in  the  two  inner  loops.  The  density  function  is  calculated  using  the 
k^,  k^,  k^,  and  k^  values  for  the  given  interaction. 

31.  The  basic  loci  contain  all  the  possible  orientations  of  k^ 
and  k^  ,  but  sometimes  this  orientation  needs  to  be  rotated  to  another 
section  of  the  grid.  Simple  rotation  procedures  can  be  used  when  the 

0  value  of  k^  is  larger  than  the  6  value  of  k^  .  Values  are  tabu¬ 
lated  in  Table  5  showing  that  loci  with  P-vectors  of  equal  magnitude  are 
equivalent  although  their  placement  in  space  is  different.  This  table 
shows  how  loci  from  interactions  with  P-vectors  of  equal  magnitude  can 
be  obtained  by  a  rotation  of  the  basic  loci.  Two  sets  of  locus  coordi¬ 
nates  are  shown — one  set  is  the  specific  iterated  coordinates  for  the 
interaction  and  the  other  set  is  the  rotation  of  the  basic  loci  coordi¬ 
nates.  Values  show  that  the  rotation  process  is  as  valid  as  iterating 
each  specific  locus.  Table  5  refers  to  rotations  where  0^  is  greater 
than  0^  ,  and  this  rotation  procedure  is  the  same  as  using  a  rotation 
matrix  on  the  basic  loci  coordinates  (in  polar  coordinates  this  reduces 
to  addition  of  a  rotation  angle). 
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Rotated  Loci 


32.  Situations  will  exist  where  the  0  value  of  is  larger 

than  the  0  value  of  k^  .  An  example  is  when  k^  has  a  0  value  of 

9°  and  k^  has  a  value  of  0°.  A  set  of  basic  loci  has  not  been  defined 

for  this  condition,  but  there  are  data  for  an  interaction  that  has  a 

^-vector  of  equal  magnitude.  Figure  8  shows  two  interactions  with  equal 

P-vectors.  The  data  for  the  ?9-locus  are  given  in  the  basic  loci  and  the 

“►  — ► 

data  for  P,  can  be  calculated  by  a  rotation  of  the  P„  data.  P 

1  2  r 

shows  the  geometrical  orientation  of  .  Table  6  shows  10  coordinates 

of  two  loci  with  P-vectors  of  the  same  magnitude.  The  third  locus  has 
been  obtained  by  using  a  rotation  procedure  on  locus  1.  The  rotation 
procedure  amounts  to  entering  the  negative  value  of  the  0^  of  locus  1 
and  adding  0^  to  the  -0^  for  the  angle  coordinate  of  the  rotated 
locus.  Comparison  of  the  locus  coordinates  of  the  rotated  locus  to 
locus  2  shows  that  they  are  the  same  locus  although  the  corresponding 
increment  numbers  may  be  different. 

33.  The  integration  arrays  are  filled  using  the  convention  that 
k^  is  always  larger  than  k^  .  To  include  the  values  for  the  transfer 
for  a  situation  where  k^  is  larger  than  k^  ,  the  program  interchanges 
the  k^  and  k^  indices  and  enters  a  negative  value  equivalent  to  the 
original  k^  -  k^  transfer  in  the  final  integral  array.  The  program 
ends  with  the  evaluation  of  dn/dt  by  using  a  summation  technique  for 
each  point  of  the  grid. 


Results 


34.  Figure  11  shows  the  contoured  ”  integration  results  for 
the  Pierson-Moskowitz  spectrum  that  Webb  obtained  using  his  program 
technique.  Figure  12  shows  the  contoured  ^  results  obtained  using 
the  integration  technique  discussed  in  the  body  of  this  report  and  the 
Pierson-Moskowitz  spectral  function.  Note  that  the  contour  plot  in  Fig¬ 
ure  12  plots  the  wave  number  (or  radial  value)  on  the  x-axis  and  the 
0-value  on  the  y-axis.  Contour  values  in  Figure  12  are  multiplied  by 
10  ^ .  Comparison  of  the  two  contours  show  similar  behavior,  and  numeri¬ 
cal  results  agree  very  closely.  Note  that  Webb's  contour  plot  uses 
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Table 


oil  0.2 J 


Figure  11.  The  nonlinear  transfer  dn/dt  as  a  function  of  wave  number 
The  contours  are  marked  in  units  of  10~^  MKS  units  (from  Webb  1978) 


OK  0  1694 


NONLINEAR  TRANSFER  FOR  PIERSON  MOSKOWITZ  SPECTRUM 


Figure  12.  Contour  plot  of  the  nonlinear  energy  transfer  for 
the  Pierson-Moskowitz  spectrum  (all  contour  values  are  multiplied 
by  10“-’)  as  a  function  of  wave-number  magnitude  and  direction 

(Tracy  1979) 


Cartesian  coordinates  for  the  wave  numbers.  The  computer  time  for  the 
Pierson-Moskowitz  spectrum  using  30  values  of  k  (0.14  m  '  to  2.44  m  S 
and  angular  increments  from  -90°  to  +90°  in  4.5°  increments  was  151  sec 
on  the  CRAY  computer. 

35.  The  Pierson-Moskowitz  spectrum  was  used  initially  +>■  results 
from  our  technique  could  be  compared  with  Webb's  results  using  the  same 
spectrum.  The  JONSWAP  researchers  list  a  general  analvtii.il  form  for  a 
spectral  density  function  (Sell  and  Hasselmann  1472).  It  contains  vari¬ 
ous  shape  parameters,  y  ,  ,  and  ,  that  can  he  varied  easily  in 

the  computer  coding.  This  spectral  density  function  can  be  written  as: 


where  cx  is  the  equilibrium  rsnge  coefficient,  f  is  the  frequency  of 

m 

the  spectral  peak,  y  is  the  peak  enhancement  parameter  and  and 

a,  are  the  left  and  right  spectral  widths  of  the  spectrum.  f  and  a 
b  m 

were  treated  as  constants  for  the  various  spectra.  This  density  function 

2 

was  used  with  a  cos  0  spreading  function.  The  Pierson-Moskowitz 
spectrum  corresponds  to  y  =  1  ,  a  =  0.07  ,  and  o,  =0.09  .  Various 

3  D 

combinations  of  shape  parameters  were  run  and  computer  times  are  similar 
to  the  run  time  for  the  Pierson-Moskowitz  spectrum. 

36.  Figures  13-17  give  the  results  of  the  calculations  of  the  net 
source  functions  for  various  spectra  of  the  type  described  by  Equation  43. 
These  results  are  given  in  terms  of  angularly  integrated  functions  since 
these  are  eacier  to  visualize  than  the  corresponding  two-dimensional 
source  functions.  Figure  13  shows  the  one-dimensional  energy  transfer 
function  for  the  Pierson-Moskowitz  spectrum  (y  =  1  ,  a  =  0.07  , 
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0.30 


frequency  for  four  JONSWAP  frequency  spectra  with  y  =  3.3  and  different  shape 

(Continued) 
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17.  The  nonlinear  energy  transfer  (S)  and  the  input  energy  E(f)  as 
function  of  frequency  for  the  Wave  Model  Spectrum  and  the  SGONO  spectrum 


a,  =  0.09)  .  The  one-dimensional  transfer  rate  is  calculated  by  S(f)  = 
S(f ,0)  d0  where  S(f,0)  is  —  and  is  graphed  in  Figure  13  and  Fig¬ 
ures  14-17  as  a  function  of  frequency.  The  spectral  energy  is  graphed 
on  the  right-hand  side  of  the  graph  and  is  also  a  function  of  frequency. 
Figure  13  can  be  compared  with  the  JONSWAP  results  for  the  Pierson- 
Moskowitz  spectrum  (Sell  and  Hasselmann  1972).  Figures  14-16  can  also 
be  compared  with  the  JONSWAP  results  (Sell  and  Hasselmann  1972).  Two 
spectra  that  do  not  correspond  to  the  JONSWAP  analytical  expression  were 
treated  in  Figure  17.  The  wave  model  spectrum  has  an  f  ^  tail  and  the 


forward  face  of  the  spectrum  is  represented  by  J exp  ]-(f  / f ) ^ 1  E,  ( 

_5  |  L  m  J  fm( 

where  f  =0.3  and  E_  a  f  .  The  SG0N0  spectrum  also  con-  ' 
m  _5  fm  51 

tains  the  f  tail  and  f  =  0.3  .  The  slope  of  the  forward  face  of 

m 

the  spectrum  was  determined  from  the  input  wind  speed  to  create  an 

f  =  0.3  .  E  is  in  M^/Hz  and  S  or  is  calculated  in  MKS  units, 

a  at 

Results  indicate  that  this  approach  to  calculation  of  the  nonlinear 
energy  transfer  agrees  with  previous  research  and  is  an  efficient 
process. 

37.  In  a  comparison  of  the  S-curves  in  Figures  13  and  14-17  a 
similar  behavior  is  observed.  This  behavior  shows  the  source-sink  effect 
of  the  nonlinear  interactions  within  the  spectrum.  The  variation  of  the 
peak  enhancement  factor,  y  ,  shows  that  as  the  spectrum  becomes  steeper 
the  negative  lobe  of  the  wave- wave  interactions  shows  a  trend  toward 
the  positive  and  in  the  very  steep  spectra  the  function  even  develops  a 
slight  positive  lobe  within  the  usual  negative  region.  In  the  steepest 
spectra  the  nonlinear  curve,  S  ,  approaches  a  sinusoidal  damping  curve. 
The  steep  spectrum  acts  like  a  single  energy  packet  at  a  certain  Af 
range  and  shows  how  the  nonlinear  reactions  attempt  to  create  a  more 
stable  spectral  form. 


biscussion 

38.  The  main  result  of  this  report  is  to  provide  an  efficient 
integration  technique  to  evaluate  Boltzmann  integrals  of  the  type  de¬ 
rived  for  energy  transfer  due  to  wave-wave  interaction  in  a  wave  spectrum. 
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In  the  past,  evaluation  of  these  integrals  required  extensive  computer 
time.  Webb's  transformation  of  the  form  of  the  integrals  simplified 
the  mathematics  but  did  not  significantly  reduce  the  computational  time 
for  their  calculation.  This  report  has  demonstrated  that  by  a  straight 
forward  transformation  of  the  integration  space,  a  large  increase  in 
computational  efficiency  is  achieved.  Webb  has  discussed  various  physi¬ 
cal  aspects  of  the  results  of  the  nonlinear  transfer  and  the  JONSWAP 
study  has  shown  these  results  in  their  experimental  data.  The  reader  is 
referred  to  these  papers  for  a  complete  discussion  of  these  results. 

39.  Specialized  numerical  techniques  and  grids  set  up  for  a 
specific  problem  are  often  only  useful  in  one  specific  problem  and  have 
no  relation  to  the  physical  problem.  The  technique  developed  in  this 
paper  sets  up  a  series  of  similar  triangles  for  a  given  0^  and  0^ 
orientation  as  shown  in  Figure  7.  This  sets  up  a  similarity  relation 
between  the  energy  transfer  of  the  various  P-vectors  shown  at  a  given 
0^  -  0^  orientation.  If  a  situation  exists  where  the  density  function 
scales,  the  nonlinear  transfer  at  a  fixed-angle  difference  could  be  eval¬ 
uated  by  a  summation.  The  approach  outlined  in  this  report  follows  the 
actual  physics  of  the  interactions  and  gives  a  geometrical  insight  into 
the  actual  process  without  obscuring  the  process  with  complex  mathemati¬ 
cal  manipulation.  These  concepts  will  be  useful  in  describing  the  com¬ 
plex  nature  of  nonlinear  wave-wave  interactions. 
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This  is  the  subroutine  used  to  calculate  the  coupling  coefficient. 
The  equation  for  the  coupling  coefficient  is  given  by  Webb  (1978).  In 
the  subroutine  below  WK1  corresponds  to  corresponds  to 
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Appendix  B:  Notation 

Dummy  variable  for  delta  function 

Group  velocity 

Coupling  coefficient 

Density  function 

Energy 

Frequency 

Frequency  of  peak  of  spectrum 

Product  of  coupling  coefficient,  0(x)  ,  and  density 
function 

Acceleration  due  to  gravity 
Geometric  product 

L  h 

Wave-number  vector  for  1  wave  number 


Smallest  value  of  r  on  (r,0)  grid 
t  h 

Action  density  of  i  wave  number 
Normal  (radial)  direction  in  k^-space 
Number  of  radial 

Interaction  vector  (k^  -  k^) 

, 1/2  . L/2  ,  .  ,  - 

k^  -  (a  magnitude) 

1*1 

Tangential  direction  in  k^-space 
Transfer  function 

Locus  equation  =  Q  +  k^2  -  ^|p  +  =  0 

Phillips  equilibrium  constant 


B1 


0 


0  (x) 


0 


A  measure  of  spectral  width  in  spectral  equation 
Peak  enhancement  parameter  in  spectral  equation 
Angular  velocity  of  iC^  wave  number 
Dirac  delta  function 

Scaling  factor  in  integral;  geometric  spacing  multiple 

0  -  k^ |  -  |£x  -  k3j)  ;  has  a  value  of  1  or  0  for 

integral 

Angle  measurement  in  polar  grid 
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